import numpy as np
import pandas as pd
from scipy import stats

import matplotlib.pyplot as plt
import japanize_matplotlib
%matplotlib inline

from tqdm import tqdm_notebook as tqdm

pd.set_option('display.max_rows', 9999)

データ整形¶

2013年9月26日の1:100の株式分割を調整したデータを用いる

df = pd.read_csv("7203_2003.csv", encoding='cp932', header=1)
for i in range(16):
    year = i+2004
    temp = pd.read_csv(f"7203_{year}.csv", encoding='cp932', header=1)
    df = pd.concat((df,temp),axis=0)
df = df.reset_index(drop=True)

df.tail()

マーケット株価としてTOPIXを採用する。

topix = pd.read_csv("topix.csv").rename(columns={'終値': 'TOPIX終値', '日付': 'TOPIX日付'})

for i in range(len(topix)):
    if (i%22 == 20) or (i%22 == 21):
        topix.drop(i, inplace=True)
topix = topix[::-1].reset_index(drop=True)
topix_copy = topix.copy()
topix.drop([3380,3382,3647,4090], axis=0, inplace=True)
topix = topix.reset_index(drop=True)

topix.head()

無リスク資産として、日本国債（10年）を採用する

bond = pd.read_csv("bond.csv")
bond.drop([3380,3382,3647,4090], axis=0, inplace=True)
bond = bond.reset_index(drop=True)

bond.tail()

date = pd.concat((df["日付"].str.replace('-',''),bond["国債日付"].str.replace('/','')),axis=1)
date.fillna(0,inplace=True)
date = date.astype(int)
for i in range(len(date)):
    if date.loc[i,"日付"]!=date.loc[i,"国債日付"]:
        print(i)

#date.loc[4142-1:4142]

#topix_copy[4087-1:4087+5]

df.shape, topix.shape, bond.shape

((4142, 7), (4142, 5), (4142, 2))

df = pd.concat((df,topix["TOPIX終値"].astype(float)),axis=1)
df = pd.concat((df,bond["金利"].astype(float)),axis=1)

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 4142 entries, 0 to 4141
Data columns (total 9 columns):
日付         4142 non-null object
始値         4142 non-null int64
高値         4142 non-null int64
安値         4142 non-null int64
終値         4142 non-null int64
出来高        4142 non-null int64
終値調整値      4142 non-null int64
TOPIX終値    4142 non-null float64
金利         4142 non-null float64
dtypes: float64(2), int64(6), object(1)
memory usage: 291.3+ KB

df["日足収益率"] = df["終値調整値"].diff()
df["TOPIX日足利益率"] = df["TOPIX終値"].diff()
df["日足収益率"][:-1] = df["日足収益率"][1:]
df["TOPIX日足利益率"][:-1] = df["TOPIX日足利益率"][1:]
df["日足収益率"] = df["日足収益率"]/df["終値調整値"]
df["TOPIX日足利益率"] = df["TOPIX日足利益率"]/df["TOPIX終値"]
df["日足収益率"][1:] = df["日足収益率"][:-1]
df["TOPIX日足利益率"][1:] = df["TOPIX日足利益率"][:-1]
df["日足収益率"][0:1]=np.nan
df["TOPIX日足利益率"][0:1]=np.nan

/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  This is separate from the ipykernel package so we can avoid doing imports until
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  after removing the cwd from sys.path.
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:7: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  import sys
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:8: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:9: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  if __name__ == '__main__':
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:10: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  # Remove the CWD from sys.path while we load stuff.

df["週足収益率"] = df["終値調整値"].diff(7)
df["TOPIX週足利益率"] = df["TOPIX終値"].diff(7)
df["週足収益率"][:-7] = df["週足収益率"][7:]
df["TOPIX週足利益率"][:-7] = df["TOPIX週足利益率"][7:]
df["週足収益率"] = df["週足収益率"]/df["終値調整値"]
df["TOPIX週足利益率"] = df["TOPIX週足利益率"]/df["TOPIX終値"]
df["週足収益率"][7:] = df["週足収益率"][:-7]
df["TOPIX週足利益率"][7:] = df["TOPIX週足利益率"][:-7]
df["週足収益率"][0:7]=np.nan
df["TOPIX週足利益率"][0:7]=np.nan

/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  This is separate from the ipykernel package so we can avoid doing imports until
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  after removing the cwd from sys.path.
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:7: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  import sys
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:8: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:9: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  if __name__ == '__main__':
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:10: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  # Remove the CWD from sys.path while we load stuff.

df["月間収益率"] = df["終値調整値"].diff(30)
df["TOPIX月間利益率"] = df["TOPIX終値"].diff(30)
df["月間収益率"][:-30] = df["月間収益率"][30:]
df["TOPIX月間利益率"][:-30] = df["TOPIX月間利益率"][30:]
df["月間収益率"] = df["月間収益率"]/df["終値調整値"]
df["TOPIX月間利益率"] = df["TOPIX月間利益率"]/df["TOPIX終値"]
df["月間収益率"][30:] = df["月間収益率"][:-30]
df["TOPIX月間利益率"][30:] = df["TOPIX月間利益率"][:-30]
df["月間収益率"][0:30]=np.nan
df["TOPIX月間利益率"][0:30]=np.nan

/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  This is separate from the ipykernel package so we can avoid doing imports until
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  after removing the cwd from sys.path.
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:7: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  import sys
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:8: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:9: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  if __name__ == '__main__':
/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:10: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  # Remove the CWD from sys.path while we load stuff.

df.tail()

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"TOPIX")
plt.plot(df["TOPIX終値"], linewidth=0.5)

[<matplotlib.lines.Line2D at 0x10a876198>]

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"日本国債・10年")
plt.plot(df["金利"], linewidth=0.5)

[<matplotlib.lines.Line2D at 0x11bb394a8>]

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"東京電力・株価")
plt.plot(df["終値調整値"], linewidth=0.5)

[<matplotlib.lines.Line2D at 0x11bd9a9b0>]

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"東京電力・株価・日足収益率")
plt.plot(df["日足収益率"][1:], linewidth=0.5)

[<matplotlib.lines.Line2D at 0x11c270278>]

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"東京電力・株価・週足収益率")
plt.plot(df["週足収益率"][7:], linewidth=0.5)

[<matplotlib.lines.Line2D at 0x11c007208>]

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"東京電力・株価・月間収益率")
plt.plot(df["月間収益率"][30:], linewidth=0.5)

[<matplotlib.lines.Line2D at 0x11c04f588>]

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"TOPIX月間超過収益率に対する東京電力の月間超過収益率")
plt.scatter(df["TOPIX月間利益率"][30:]-df["金利"][30:],df["月間収益率"][30:]-df["金利"][30:], s=5)
plt.xlabel('TOPIX月間収益率')
plt.ylabel('月間収益率')
plt.grid()
plt.show()

x = df["月間収益率"][30:]-df["金利"][30:]
length = len(x)

x.shape

(4112,)

x.mean()

-0.8743035933088424

x.std()

0.6319247506951736

MCMC¶

epsilon = 0.01
# update width
trial = 25000
# trial: how many times do mcmc

theta = [0.9,1.1,0.001,0.005,1,1,0,0]
theta[6] = np.exp(theta[4])/(1+np.exp(theta[4]))
theta[7] = np.exp(theta[5])/(1+np.exp(theta[5]))
# theta[0] = beta1
# theta[1] = beta2
# theta[2] = sigma1
# theta[3] = sigma2
# theta[4] = p0
# theta[5] = q0
# theta[6] = p
# theta[7] = q

num_val = len(theta)

対数尤度を計算する関数

excess_return = df["TOPIX月間利益率"][30:]-df["金利"][30:]

def cal_logL(x, theta):
    likelihood1 = stats.norm.pdf(x=x, loc=theta[0]*excess_return, scale=np.sqrt(theta[2]))
    likelihood2 = stats.norm.pdf(x=x, loc=theta[1]*excess_return, scale=np.sqrt(theta[3]))
    
    prior1=(1-theta[7])/(2-theta[6]-theta[7])
    prior2=(1-theta[6])/(2-theta[6]-theta[7])
    
    logL = 0
    for i in range(length):
        temp1 = likelihood1[i]*prior1
        temp2 = likelihood2[i]*prior2
        
        if (temp1 == 0) & (temp2 == 0):
            temp1 +=  10**(-100)
            temp2 +=  10**(-100)

        posterior1 = temp1/(temp1+temp2)
        posterior2 = temp2/(temp1+temp2)

        logL += np.log(temp1+temp2)
        
        prior1 = theta[6]*posterior1 + (1-theta[7])*posterior2
        prior2 = (1-theta[6])*posterior1 + theta[7]*posterior2
        
    return logL

各レジームである確率の配列を返す関数

def prob_regime(x, theta):
    likelihood1 = stats.norm.pdf(x=x, loc=theta[0], scale=np.sqrt(theta[2]))
    likelihood2 = stats.norm.pdf(x=x, loc=theta[1], scale=np.sqrt(theta[3]))
    
    prior1=(1-theta[7])/(2-theta[6]-theta[7])
    prior2=(1-theta[6])/(2-theta[6]-theta[7])
    
    prob_list1 = []
    prob_list2 = []
    for i in range(length):
        temp1 = likelihood1[i]*prior1
        temp2 = likelihood2[i]*prior2
        
        if (temp1 == 0) & (temp2 == 0):
            temp1 +=  10**(-100)
            temp2 +=  10**(-100)

        prob_list1.append(temp1/(temp1+temp2))
        prob_list2.append(temp2/(temp1+temp2))
        
        posterior1 = temp1/(temp1+temp2)
        posterior2 = temp2/(temp1+temp2)
        
        prior1 = theta[6]*posterior1 + (1-theta[7])*posterior2
        prior2 = (1-theta[6])*posterior1 + theta[7]*posterior2
    
    return prob_list1, prob_list2

$\theta$を更新する関数。
レジーム1よりレジーム2のボラティリティが大きいことを条件としている。

def create_next_theta(theta, epsilon):
    new_theta = theta.copy()
    new_theta[0:6] += (2*np.random.rand(6)-1)*epsilon
    new_theta[2:4] = np.exp(np.log(theta[2:4]) + (2*np.random.rand(2)-1)*epsilon*0.1)
    
    if new_theta[0]>new_theta[1]:
        new_theta[0] = theta[0]
        new_theta[1] = theta[1]

    new_theta[6] = np.exp(new_theta[4])/(1+np.exp(new_theta[4]))
    new_theta[7] = np.exp(new_theta[5])/(1+np.exp(new_theta[5]))
    
    return list(new_theta)

MCMCを実行する関数

def mcmc(x, theta, trial):
    # x: time_series
    # theta: parameter
    # trial: how many times do mcmc
    
    theta_list = []
    logL_list = []
    theta_list.append(theta)
    
    for i in tqdm(range(trial)):
        next_theta = create_next_theta(theta, epsilon)
        
        logL = cal_logL(x, theta)
        next_logL = cal_logL(x, next_theta)
        
        ratio = np.exp(next_logL-logL)
        logL_list.append(logL)

        if ratio > 1:
            theta = next_theta

        elif ratio > np.random.rand():
            theta = next_theta

        theta_list.append(theta)
        
        #if np.abs(np.log(ratio)) < 10**(-8):
            #return theta_list
        
        #if (i <120) & (i%4==0):
        if (False):
            prob1, prob2 = prob_regime(x, theta)
            regime = np.array(prob2)>=0.5
            regime = regime.astype(int)

            fig = plt.figure(figsize=(10,5),dpi=200)
            ax1 = fig.add_subplot(111)
            ln1=ax1.plot(range(length), x, linewidth=0.7, color="black", label="data")

            ax2 = ax1.twinx()
            ln2=ax2.bar(range(length), 1-np.array(regime), width=1.0, color="blue", alpha=0.4, label=r"$P(I(t) = 1)$")
            ln2=ax2.bar(range(length), regime, width=1.0, color="red", alpha=0.4, label=r"$P(I(t) = 2)$")

            h1, l1 = ax1.get_legend_handles_labels()
            h2, l2 = ax2.get_legend_handles_labels()
            ax1.legend(h1+h2, l1+l2, bbox_to_anchor=(1.05, 1), loc='upper left', edgecolor="gray", fancybox=False)

            ax1.set_xlabel('t')
            ax1.set_ylabel(r'data')
            ax2.set_ylabel(r'$P(I(t) = 2)$')
            ax1.set_title(f"Estimated Regime, epsilon={epsilon}, trial={i}/{trial}")
            plt.savefig(f'regime_estimating_{i}.png')
            
        if i%100==0:
            print(logL)
    
    return theta_list, logL_list

実行¶

theta_list, logL_list = mcmc(x, theta, trial)

2255.826130746799

/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:16: RuntimeWarning: overflow encountered in exp
  app.launch_new_instance()

6722.651502517366
6742.561729442786
6757.4465031138125
6788.5121950481935
6815.652893017021
6838.561493789909
6845.173400920939
6863.888546392502
6877.231776113678
6894.669135018193
6913.3486517095735
6929.515925873457
6932.474675405062
6953.324375939682
6972.583859922756
6990.038093795388
7005.476537174056
7012.1710357810025
7020.001448235798
7029.089915945403
7041.862431272022
7048.907478764396
7054.371921138474
7062.71119057893
7068.172619188613
7070.389869551357
7075.847652390752
7078.734079124669
7080.4578177028825
7083.091130095027
7084.770656836493
7093.290574131
7094.452076931246
7100.801504488748
7107.436202951625
7107.654535514021
7111.113770269181
7114.065877915383
7119.378653528787
7123.358209435673
7127.90225902587
7126.5007576402
7132.281164140301
7138.308291626326
7139.580545250712
7139.182744572483
7141.409878929371
7144.31902633542
7146.68865102671
7147.057155545789
7147.748768400299
7147.567152982717
7147.894970222863
7152.359106431969
7153.119213995057
7155.24163372912
7156.99159060253
7157.988415526203
7158.374178984363
7159.0526676010295
7160.879392692513
7162.373890138914
7163.052072472195
7162.039518149662
7160.92518187674
7163.009201341686
7163.598766419299
7164.746871667794
7165.605007406617
7165.976126396313
7166.700863352339
7165.3519834119825
7166.48959696817
7167.64101646207
7167.355041481791
7169.668173550249
7169.631180060298
7171.969906701699
7173.247720584662
7171.726202985625
7172.451780514191
7173.380464755888
7173.070371634205
7174.23722347523
7176.930520527941
7177.761001213441
7176.162500762713
7176.776205324267
7177.763483088963
7178.4846471966675
7180.030660651551
7179.260992986536
7180.535848224759
7181.789854578718
7184.295356567734
7185.028161350145
7183.34468535788
7182.486587560096
7182.587964588489
7181.83859782041
7182.671092609848
7182.528069530775
7185.314372855485
7186.716436752039
7186.2479442896065
7187.290633517316
7189.494820477167
7190.477971342918
7189.349263744632
7189.805795372553
7188.035279615646
7187.884139769337
7188.89983525508
7188.746100961056
7188.159486840969
7186.963571731484
7191.422549621078
7191.432589794269
7193.292286295713
7193.782233922268
7195.630003773168
7196.9003966157825
7199.533593662323
7198.399121865576
7199.570848591694
7199.632248775376
7199.114027497195
7199.268595446606
7200.136158549556
7198.466157351381
7197.705043985156
7198.887469982997
7198.320346793627
7197.126982538791
7196.325756780242
7197.678023863074
7197.340094658124
7199.286777786894
7199.337456444561
7198.355102904891
7199.294197764726
7197.598184440031
7194.429008378105
7198.326208490215
7197.436008227187
7199.368497388582
7200.174329565148
7202.0580367866305
7203.843779709368
7206.854963555146
7203.866299048918
7206.219410273423
7206.775936124377
7208.603106040562
7208.761659443578
7209.7346862302165
7209.299530642327
7209.479192439506
7213.6054470181225
7215.476494767593
7216.258342240586
7215.126821155817
7216.786373747356
7216.591907797151
7217.926142296173
7217.8191818814585
7219.779672993653
7219.880309210739
7220.618244157974
7220.877320836695
7220.139676347293
7219.275207198224
7219.454572172441
7219.327682674991
7219.537527009822
7218.5293569776495
7218.930461069524
7219.8096865410935
7219.559740308576
7219.323590184995
7220.735974380344
7221.192107589889
7217.742691336136
7217.719565311187
7219.516643219143
7218.382192398013
7218.289285591092
7218.894114170123
7219.814950888878
7224.6519208491645
7223.552307345189
7224.238164415689
7224.213388135138
7225.564663076441
7225.209611064689
7225.2514168086345
7226.344288780983
7226.493638186632
7229.382704564211
7228.628138977993
7228.767146543674
7229.688645310213
7230.165955833863
7230.698517585316
7232.460133891387
7232.015220823591
7229.780324264314
7230.85482426318
7232.753636918729
7234.004599548687
7234.611459928468
7233.339200284027
7232.808801312365
7233.527039902251
7233.164996241022
7233.664943141734
7234.312348704496
7235.040431236862
7235.417947691474
7234.387165291439
7235.536034369914
7236.401142144173
7235.495583971092
7236.179059013223
7237.40792847301
7237.709429807006
7237.7461910045
7236.63219545131
7235.941643306059
7237.086326670372
7237.81478899965
7239.208500351918
7239.790975424525
7238.797812553145
7237.694829938112
7239.331601679291
7240.457421721558
7241.381423489751
7237.438649505965
7238.817362159183
7237.339559088966
7238.592836103577
7240.770764211652
7240.297601281493
7243.142063836352
7244.028880009606
7244.72120579173
7245.091506881953
7245.600066903025

theta_list = np.array(theta_list)

prob1, prob2 = prob_regime(x, theta_list[-1])

theta_list[-1]
# theta[0] = beta1
# theta[1] = beta2
# theta[2] = sigma1
# theta[3] = sigma2
# theta[4] = p0
# theta[5] = q0
# theta[6] = p
# theta[7] = q

array([1.00156068e+00, 1.00183292e+00, 6.38503075e-04, 5.32704755e-03,
       3.62403199e+00, 3.28519943e+00, 9.74018157e-01, 9.63917559e-01])

prob = pd.DataFrame(np.round(prob1,3))
prob.index = prob.index + 30
pd.concat((df,prob),axis=1)

可視化¶

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"月間超過収益率")
plt.plot(x, linewidth=0.5, label="トヨタ自動車")
plt.plot(excess_return, linewidth=0.5, label="TOPIX")
plt.legend(edgecolor="gray", fancybox=False)
plt.xlabel('t')
plt.ylabel(r'超過収益率')
plt.savefig("excess_return_series.png")

plt.figure(figsize=(10,5),dpi=200)
plt.title(f"TOPIX月間超過収益率に対するトヨタ自動車の月間超過収益率")
plt.scatter(df["TOPIX月間利益率"][30:]-df["金利"][30:],df["月間収益率"][30:]-df["金利"][30:], s=5)
plt.xlabel('TOPIX月間収益率')
plt.ylabel('月間収益率')
plt.grid()
plt.savefig("excess_return_scatter.png")
plt.show()

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ln1=ax1.plot(range(length), x, linewidth = 1.0, label="return")

ax2 = ax1.twinx()
ln2=ax2.bar(range(length), prob2, width=1.0, color="red", alpha=0.4, label=r"$P(I(t) = 2)$")

h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, bbox_to_anchor=(1.05, 1), loc='upper left')

ax1.set_xlabel('t')
ax1.set_ylabel(r'超過収益率')
ax2.set_ylabel(r'$P(I(t) = 2)$')
ax1.set_title(f"トヨタ自動車・高ベータレジームである確率, epsilon={epsilon}, trial={trial}")

plt.subplots_adjust(left = 0.1, right = 0.8)
plt.savefig("probability_regime2.png")
plt.show()

regime = np.array(prob2)>=0.5
regime = regime.astype(int)

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ln1=ax1.plot(range(length), x, linewidth=0.7, color="black", label="return")

ax2 = ax1.twinx()
ln2=ax2.bar(range(length), 1-np.array(regime), width=1.0, color="blue", alpha=0.4, label=r"$P(I(t) = 1)$")
ln2=ax2.bar(range(length), regime, width=1.0, color="red", alpha=0.4, label=r"$P(I(t) = 2)$")

h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, bbox_to_anchor=(1.05, 1), loc='upper left', edgecolor="gray", fancybox=False)

ax1.set_xlabel('t')
ax1.set_ylabel(r'超過収益率')
ax1.set_title(f"トヨタ自動車・予測レジーム, epsilon={epsilon}, trial={trial}")

plt.subplots_adjust(left = 0.1, right = 0.8)
plt.savefig("estimated_regime_with_return.png")
plt.show()

regime = np.array(prob2)>=0.5
regime = regime.astype(int)

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ln1=ax1.plot(range(length), df["終値調整値"][30:], linewidth=0.7, color="black", label="stock price")

ax2 = ax1.twinx()
ln2=ax2.bar(range(length), 1-np.array(regime), width=1.0, color="blue", alpha=0.4, label=r"$P(I(t) = 1)$")
ln2=ax2.bar(range(length), regime, width=1.0, color="red", alpha=0.4, label=r"$P(I(t) = 2)$")

h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, bbox_to_anchor=(1.05, 1), loc='upper left', edgecolor="gray", fancybox=False)

ax1.set_xlabel('t')
ax1.set_ylabel(r'株価')
ax1.set_title(f"トヨタ自動車・予測レジーム, epsilon={epsilon}, trial={trial}")

plt.subplots_adjust(left = 0.1, right = 0.8)
plt.savefig("estimated_regime_with_price.png")
plt.show()

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ax1.plot(theta_list[:,6], linewidth=0.5, color="navy", label=r"$P(I(t) = 1 | I(t-1) = 1)$")
ax1.plot(theta_list[:,7], linewidth=0.5, color="darkred", label=r"$P(I(t) = 2 | I(t-1) = 2)$")
ax1.legend(edgecolor="gray", fancybox=False)
ax1.set_title(f"トヨタ自動車・同一レジームの滞留確率の更新過程, epsilon={epsilon}, trial={trial}")
ax1.set_xlabel('trial')
ax1.set_ylabel(r'確率')

plt.savefig("updating_probability.png")
plt.show()

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ax1.plot(theta_list[:,0], linewidth=0.5, color="navy", label=r"$\beta_1$")
ax1.plot(theta_list[:,1], linewidth=0.5, color="darkred", label=r"$\beta_2$")
ax1.legend(edgecolor="gray", fancybox=False)
ax1.set_title(f"トヨタ自動車・ベータの更新過程, epsilon={epsilon}, trial={trial}")
ax1.set_xlabel('trial')
ax1.set_ylabel(r'$\beta$')

plt.savefig("updating_mu.png")
plt.show()

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ax1.plot(theta_list[:,2], linewidth=0.5, color="navy", label=r"$\sigma_1$")
ax1.plot(theta_list[:,3], linewidth=0.5, color="darkred", label=r"$\sigma_2$")
ax1.legend(edgecolor="gray", fancybox=False)
ax1.set_title(f"トヨタ自動車・$\sigma$の更新過程, epsilon={epsilon}, trial={trial}")
ax1.set_xlabel('trial')
ax1.set_ylabel(r'$\sigma$')

plt.savefig("updating_sigma.png")
plt.show()

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ax1.plot(logL_list, linewidth=0.5)
ax1.set_title(f"対数尤度の推移")
plt.savefig("likelihood.png")
plt.show()

	日付	始値	高値	安値	終値	出来高	終値調整値
4137	2019-11-26	7788	7820	7707	7708	4567100	7708
4138	2019-11-27	7706	7728	7685	7690	2744100	7690
4139	2019-11-28	7740	7757	7688	7724	3226900	7724
4140	2019-11-29	7740	7740	7638	7638	3762600	7638
4141	2019-12-02	7701	7761	7682	7738	3913300	7738

	TOPIX日付	始値	高値	安値	TOPIX終値
0	2003/01/06	851.8	862.19	851.8	860.07
1	2003/01/07	868.32	870.69	853.93	853.93
2	2003/01/08	850.05	850.05	839.55	839.55
3	2003/01/09	832.69	839.33	829.69	839.33
4	2003/01/10	843.31	843.83	830.92	837.7

	国債日付	金利
4137	2019/11/26	-0.092
4138	2019/11/27	-0.106
4139	2019/11/28	-0.086
4140	2019/11/29	-0.075
4141	2019/12/02	-0.046

	日付	始値	高値	安値	終値	出来高	終値調整値	TOPIX終値	金利	日足収益率	TOPIX日足利益率	週足収益率	TOPIX週足利益率	月間収益率	TOPIX月間利益率
4137	2019-11-26	7788	7820	7707	7708	4567100	7708	1705.71	-0.092	-0.002717	0.001615	-0.015078	0.005328	0.084412	0.078594
4138	2019-11-27	7706	7728	7685	7690	2744100	7690	1710.98	-0.106	-0.002335	0.003090	-0.018632	0.006033	0.057917	0.072533
4139	2019-11-28	7740	7757	7688	7724	3226900	7724	1708.06	-0.086	0.004421	-0.001707	-0.003869	0.006678	0.045904	0.054228
4140	2019-11-29	7740	7740	7638	7638	3762600	7638	1699.36	-0.075	-0.011134	-0.005093	-0.006504	0.004878	0.024547	0.041587
4141	2019-12-02	7701	7761	7682	7738	3913300	7738	1714.49	-0.046	0.013092	0.008903	0.002851	0.014863	0.048083	0.055616

	日付	始値	高値	安値	終値	出来高	終値調整値	TOPIX終値	金利	日足収益率	TOPIX日足利益率	週足収益率	TOPIX週足利益率	月間収益率	TOPIX月間利益率	0
0	2003-01-06	3250	3270	3240	3250	3531400	3250	860.07	0.911	NaN	NaN	NaN	NaN	NaN	NaN	NaN
1	2003-01-07	3290	3300	3220	3240	5907800	3240	853.93	0.906	-0.003077	-0.007139	NaN	NaN	NaN	NaN	NaN
2	2003-01-08	3230	3230	3160	3170	3242100	3170	839.55	0.892	-0.021605	-0.016840	NaN	NaN	NaN	NaN	NaN
3	2003-01-09	3150	3180	3130	3150	3947300	3150	839.33	0.852	-0.006309	-0.000262	NaN	NaN	NaN	NaN	NaN
4	2003-01-10	3160	3170	3100	3140	4010200	3140	837.70	0.870	-0.003175	-0.001942	NaN	NaN	NaN	NaN	NaN
5	2003-01-14	3130	3140	3080	3100	4175800	3100	845.90	0.872	-0.012739	0.009789	NaN	NaN	NaN	NaN	NaN
6	2003-01-15	3120	3120	3090	3120	4586000	3120	855.66	0.849	0.006452	0.011538	NaN	NaN	NaN	NaN	NaN
7	2003-01-16	3100	3110	3070	3090	4417500	3090	856.30	0.840	-0.009615	0.000748	-0.049231	-0.004383	NaN	NaN	NaN
8	2003-01-17	3050	3090	3050	3060	4945200	3060	859.25	0.834	-0.009709	0.003445	-0.055556	0.006230	NaN	NaN	NaN
9	2003-01-20	3030	3040	2990	3010	5897700	3010	853.90	0.825	-0.016340	-0.006226	-0.050473	0.017092	NaN	NaN	NaN
10	2003-01-21	3010	3070	3010	3050	6200100	3050	863.64	0.825	0.013289	0.011406	-0.031746	0.028964	NaN	NaN	NaN
11	2003-01-22	3020	3030	2990	3000	6093300	3000	854.06	0.811	-0.016393	-0.011093	-0.044586	0.019530	NaN	NaN	NaN
12	2003-01-23	3010	3040	2995	3040	5011800	3040	865.43	0.806	0.013333	0.013313	-0.019355	0.023088	NaN	NaN	NaN
13	2003-01-24	3010	3030	2980	2985	6135400	2985	861.36	0.812	-0.018092	-0.004703	-0.043269	0.006662	NaN	NaN	NaN
14	2003-01-27	2980	2990	2960	2975	3969700	2975	850.03	0.802	-0.003350	-0.013154	-0.037217	-0.007322	NaN	NaN	NaN
15	2003-01-28	2965	2970	2940	2945	4985300	2945	841.12	0.791	-0.010084	-0.010482	-0.037582	-0.021100	NaN	NaN	NaN
16	2003-01-29	2950	2950	2880	2895	5550500	2895	825.72	0.781	-0.016978	-0.018309	-0.038206	-0.033002	NaN	NaN	NaN
17	2003-01-30	2855	2880	2840	2865	6884400	2865	825.16	0.759	-0.010363	-0.000678	-0.060656	-0.044556	NaN	NaN	NaN
18	2003-01-31	2870	2870	2835	2855	6536500	2855	821.18	0.810	-0.003490	-0.004823	-0.048333	-0.038498	NaN	NaN	NaN
19	2003-02-03	2860	2960	2855	2955	8112700	2955	837.70	0.817	0.035026	0.020117	-0.027961	-0.032042	NaN	NaN	NaN
20	2003-02-04	2995	3000	2975	2975	6465300	2975	844.26	0.856	0.006768	0.007831	-0.003350	-0.019852	NaN	NaN	NaN
21	2003-02-05	2955	3010	2950	2965	3918200	2965	845.29	0.861	-0.003361	0.001220	-0.003361	-0.005576	NaN	NaN	NaN
22	2003-02-06	3010	3010	2965	2965	5010200	2965	839.95	0.857	0.000000	-0.006317	0.006791	-0.001391	NaN	NaN	NaN
23	2003-02-07	2970	2985	2945	2975	3276300	2975	839.11	0.839	0.003373	-0.001000	0.027634	0.016216	NaN	NaN	NaN
24	2003-02-10	2965	2975	2945	2965	2242500	2965	842.69	0.821	-0.003361	0.004266	0.034904	0.021244	NaN	NaN	NaN
25	2003-02-12	3000	3000	2970	2990	4564300	2990	857.23	0.841	0.008432	0.017254	0.047285	0.043900	NaN	NaN	NaN
26	2003-02-13	2970	2980	2930	2965	3174000	2965	853.01	0.838	-0.008361	-0.004923	0.003384	0.018276	NaN	NaN	NaN
27	2003-02-14	2935	2990	2935	2975	5047600	2975	858.05	0.839	0.003373	0.005908	0.000000	0.016334	NaN	NaN	NaN
28	2003-02-17	2990	3010	2960	2980	5143600	2980	861.70	0.870	0.001681	0.004254	0.005059	0.019413	NaN	NaN	NaN
29	2003-02-18	2965	2980	2905	2910	4776500	2910	856.70	0.872	-0.023490	-0.005802	-0.018550	0.019942	NaN	NaN	NaN
30	2003-02-19	2950	2960	2910	2935	3629600	2935	853.59	0.863	0.008591	-0.003630	-0.013445	0.017256	-0.096923	-0.007534	0.0
31	2003-02-20	2930	2950	2920	2950	2365600	2950	849.35	0.857	0.005111	-0.004967	-0.005059	0.007903	-0.089506	-0.005363	0.0
32	2003-02-21	2915	2930	2880	2885	4648900	2885	840.11	0.833	-0.022034	-0.010879	-0.035117	-0.019971	-0.089905	0.000667	0.0
33	2003-02-24	2895	2905	2870	2875	3868400	2875	838.96	0.807	-0.003466	-0.001369	-0.030354	-0.016471	-0.087302	-0.000441	0.0
34	2003-02-25	2850	2860	2795	2805	5971600	2805	819.18	0.776	-0.024348	-0.023577	-0.057143	-0.045300	-0.106688	-0.022108	0.0
35	2003-02-26	2805	2830	2805	2825	4341000	2825	818.38	0.791	0.007130	-0.000977	-0.052013	-0.050273	-0.088710	-0.032533	0.0
36	2003-02-27	2820	2825	2790	2805	3447700	2805	819.20	0.791	-0.007080	0.001002	-0.036082	-0.043773	-0.100962	-0.042610	0.0
37	2003-02-28	2825	2835	2760	2775	4770700	2775	818.73	0.786	-0.010695	-0.000574	-0.054514	-0.040839	-0.101942	-0.043875	0.0
38	2003-03-03	2815	2860	2810	2850	5987800	2850	825.63	0.787	0.027027	0.008428	-0.033898	-0.027927	-0.068627	-0.039127	0.0
39	2003-03-04	2850	2875	2850	2865	5173200	2865	828.05	0.748	0.005263	0.002931	-0.006932	-0.014355	-0.048173	-0.030273	0.0